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ABSTRACT 


There  is  a  general  consensus  that  3D  reference  models  can  be  used  to  isolate  effects  of  wave  propagation  and  thus 
help  in  improving  characterization  of  seismic  sources.  Advances  in  computation  and  numerical  method  have  made  it 
possible  to  capture  increasingly  broadband,  full  wave  generation  and  propagation  in  3D  earth  models.  The  main 
objectives  of  this  project  are  to  construct  hierarchical,  multi-grid  (resolution)  joint  P  and  S  velocity  models  and  the 
corresponding  finite-difference  strain  Green’s  tensors  (FDSGT)  for  rapid  seismic  moment  determination.  In  the  first 
year  of  a  three-year  project,  we  collected  and  processed  up  to  21  years  continuous  seismic  data  from  broadband 
seismic  stations  in  the  eastern  hemisphere  (latitude  55°S-55°N;  longitude  30°W-157°E).  Empirical  Green’s 
functions  derived  from  ambient  seismic  noise  show  clear  Rayleigh  waves  over  a  broad  frequency  range.  We 
measured  over  ~1 0,000  phase  delays  in  6  period  bands  ranging  from  50  s  to  600  s,  which  provide  new  constraints  on 
the  entire  upper  mantle,  including  the  upper  mantle  transition  zone.  Several  iterations  of  full-wave  tomographic 
inversion  have  been  carried  out  to  obtain  a  joint  P  and  S  velocity  model.  Preliminary  results  clearly  show 
high-velocity  anomalies  associated  with  plate  subduction  beneath  Indonesia,  southern  Tibet,  Iran,  and  the  Hellenic 
arc.  The  African  cratons  and  the  west  Australian  craton  appear  as  prominent  high-velocity  features.  The  eastern 
African  Rift  and  Afar  have  deep-rooted  low- velocity  anomalies.  Future  work  includes  integration  of  earthquake 
arrivals  and  ambient  noise  data. 
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OBJECTIVES 


The  main  objectives  of  this  project  are  to  construct  hierarchical,  multi-grid  (resolution)  joint  P  and  S  velocity  models 
of  the  Middle  East  and  East  Eurasia  and  the  corresponding  FDSGTs  for  rapid  earthquake/explosion  moment 
determination. 

RESEARCH  ACCOMPLISHED 


Much  of  our  effort  during  the  first  year  of  this  three-year  project  have  been  on  the  development  of  a  full- wave 
tomographic  method  using  ambient  seismic  noise  and  its  application  to  obtain  a  joint  P  and  S  velocity  model  of  the 
eastern  hemisphere  (latitude  55°S-55°N,  longitude  30°W-157°E).  In  the  following,  we  begin  with  a  brief 
introduction  to  the  project  and  the  approach  we  are  taking.  We  report  the  preliminary  results  of  full-wave  ambient 
noise  tomography  of  the  eastern  hemisphere  and  conclude  with  plans  for  remaining  work. 

Introduction 

There  is  a  general  consensus  that  the  use  of  3D  models  improves  seismic  characterization  (Antoun  et  al.,  2008). 
Advances  in  computation  and  numerical  methods  have  made  it  possible  to  capture  increasingly  broadband,  full  wave 
generation  and  propagation  in  3D  Earth  models.  But  it  is  still  impractical  to  routinely  simulate  forward  wave 
propagation  from  the  source  in  3D  regional  models  in  real  time.  One  way  to  overcome  this  computational  challenge 
is  to  take  advantage  of  an  important  property  in  seismic  wave  propagation,  the  source-receiver  reciprocity.  The  use 
of  the  reciprocity  property  and  Green’s  functions  can  drastically  reduce  computational  costs  when  the  sources 
outnumber  the  stations  (e.g.,  Eisner  and  Clayton,  2001;  Graves  and  Wald,  2001).  More  importantly,  with  3D 
Green’s  functions  pre-calculated  and  saved  in  a  database,  they  can  be  extracted  quickly  from  the  database  after  a 
seismic  event  and  used  to  determine  seismic  moments  (Zhao  et  al.,  2006;  Shen  et  al.,  201 1). 

Figure  1  is  a  schematic  illustration  of  the  unified  tomography  and  source  moment  inversion  based  on  3D  Strain 
Green  Tensor  (SGT)  databases.  Station  SGTs  are  constructed  from  3D  reference  models  by  finite-difference 
simulation  of  the  responses  to  orthogonal  unit  impulsive  point  forces  acting  at  stations.  By  extracting  SGTs  in  a 
small  volume  surrounding  the  source  reference  location,  we  can  invert  for  source  moment  tensors  and  location  in  a 
global  optimization  scheme.  Travel  time  and  amplitude  anomalies  are  measured  from  the  observed  and  synthetic 
waveforms  at  stations.  The  forward  wave  field  from  the  source  and  SGTs  are  used  to  calculate  finite-frequency 
sensitivities  of  the  measurements  to  perturbations  in  Vp,  Vs  (Zhao  et  al.,  2005;  Zhang  et  al.,  2007;  Zhang  and  Shen, 
2008;  Shen  et  al.,  2008).  Together  with  the  sensitivities  to  source  parameters,  the  measurements  and  sensitivity 
kernels  are  used  to  invert  for  the  Earth’s  structure  (e.g.,  Chen  et  al.,  2007;  Shen  and  Zhang,  2010).  The  inversion 
results  and  additional  constraints  (e.g.,  receiver  function  solutions)  are  added  to  the  3D  model.  This  process  can  be 
repeated  to  progressively  improve  the  solutions.  For  full- wave  tomography  using  ambient  noise,  the  source  moment 
inversion  is  omitted  and  the  synthetic  waveforms  at  stations  and  sensitivity  kernels  are  all  calculated  from  SGTs. 


3D  structural  models 


Figure  1.  Flow  chart  of  the  unified  tomography  and  moment  tensor  inversions  based  on  3D  FDSGT 
databases. 
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One  of  the  main  challenges  in  full-wave  tomography  is  the  intense  computation  associated  with  simulating  wave 
propagation  in  3D  models.  To  reduce  computational  cost,  which  often  scales  linearly  with  the  number  of  discrete 
nodes  in  numerical  methods,  we  adapt  a  multigrid/multilevel  method  (Briggs,  1987)  to  solve  the  wave  propagation 
and  inversion  problem  using  a  hierarchy  of  discretization  (Figure  2).  We  use  three  levels  of  finite-difference 
computational  grids  at  the  hemispherical,  regional-,  and  local-scales  to  progressively  construct  a  hierarchical, 
multi-resolution  model.  Data  analysis  will  correspondingly  evolve  from  long  period  waves  for  teleseismic 
observations  to  broadband  waves  for  local  and  regional  observations.  Two  areas  are  selected  for  detailed  studies: 
(1)  the  Middle  East,  and  (2)  the  Tibetan  plateau  and  surrounding  areas.  Each  focus  area  will  have  a  self-consistent 
P  and  S  velocity  model  from  full-wave  tomography  and  the  FDSGT  databases  for  seismic  stations  at  local  and 
teleseismic  distances. 


longitude 


Figure  2.  A  hierarchical,  multi-grid  (resolution)  approach  to  develop  a  joint  P  and  S  velocity  model  in 

Eurasia,  Africa  and  Australia.  The  three  levels  of  grid  represent  the  finite-difference  grids  used  in 
simulations  of  wave  propagation  in  3D  models  at  continental,  regional,  and  local  scales.  The  finest 
grid  has  a  grid  spacing  of  500-1000  m.  For  legibility,  the  grids  shown  have  been  down-sampled  to 
approximately  every  25  grids.  Triangles  mark  the  permanent  (red)  and  temporary  (blue)  stations 
that  are  available  or  will  be  available  during  the  project  at  the  IRIS  DMC.  Also  plotted  are  the 
broadband  Japanese  seismic  network  (F-net,  green  triangles),  Iranian  National  Seismic  Network 
(INSN,  green  triangles),  and  Iranian  Seismic  Telemetry  Network  (IRSC,  green  crosses). 

Full-Wave  Ambient  Noise  Tomography  of  the  Eastern  Hemisphere 

Reference  Model  and  Wave  Propagation  Simulation 

The  initial  3D  reference  model  is  composed  of  CRUST  2.0  (Bassin  et  al.,  2000)  for  the  crust  and  AK135 
(Kennett  et  al.,  1995)  for  the  mantle.  The  choice  of  an  average  Earth  model  for  the  mantle  instead  of  one  of  the 
current  3D  models  (e.g.,  Ritsema  and  van  Heijst,  2000;  Ritzwoller  et  al.,  2002;  Antolik  et  al.,  2003)  is  a  tradeoff 
between  finding  a  good  starting  model  and  the  needs  to  explore  a  wide  enough  model  space  and  validate  the  final 
solution. 

We  use  a  non-staggered  grid  to  simulate  wave  propagation  in  a  spherical  coordinate  (Figure  3,  Zhang  et  al.,  2011). 
We  implement  the  DRP/opt  MacCormack  scheme,  which  alternately  uses  forward  and  backward  finite-difference 
operators,  in  the  calculation.  The  solution  is  optimized  for  both  dispersion  error  and  dissipation  error,  and  has  a 
4th-order  dispersion  accuracy.  In  order  to  increase  the  efficiency  of  the  finite-difference  algorithm,  we  use  a  grid 
with  non-uniform  grid  spacing  to  discretize  the  computational  domain.  The  grid  varies  continuously  with  smaller 
spacing  in  low  velocity  layers  and  with  larger  spacing  otherwise.  The  computational  domain  is  surrounded  by 
complex-frequency  shifted  Perfect  Matched  Layers  (PMLs)  implemented  through  auxiliary  differential  equations 
(Zhang  and  Shen,  2010). 
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Figure  3.  Geometry  of  non-uniform,  non-staggered  finite-difference  grid  in  a  spherical  coordinate  system. 


The  computational  domain  covers  the  eastern  hemisphere  (latitude  55°S-55°N;  longitude  30°W-157°E)  with  a 
lateral  grid  spacing  of  0.2°,  and  a  non-uniform  grid  from  the  surface  to  -2100  km  depth  (3  km  near  surface,  and 
gradually  increasing  to  -40  km  near  the  bottom  of  the  model).  There  are  12  PMLs  surrounding  the  model 
boundaries  except  the  free  surface.  Figure  4  shows  a  snapshot  of  simulated  seismic  response  to  a  vertical  force  at 
station  BJI  (Beijing). 


Figure  4.  A  snapshot  of  the  seismic  response  at  the  surface  (vertical  ground-motion  velocity)  to  a  vertical 
force  at  station  BJT  (Beijing). 
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Empirical  Green’s  Functions 

We  collected  continuous  data  from  broadband  seismic  stations  in  the  eastern  hemisphere  during  1990  to  2010.  We 
selected  184  stations  having  at  least  300  days  of  continuous  records  for  further  analysis.  We  removed  instrument 
response  from  the  data  and  normalized  waveforms  by  the  Time  Frequency  Normalization  method  (Ekstrom  et  al., 
2009).  We  removed  time  windows  with  large  (Ms>5.5)  earthquakes.  Daily  records  from  station  pairs  were 
cross-correlated  and  stacked  to  obtain  the  Empirical  Green’s  Functions  (EGF).  Figure  5  shows  an  example  of  EGFs 
from  cross  correlation  of  vertical- vertical  component.  Our  results  confirm  that  Rayleigh  waves  can  be  clearly 
extracted  from  seismic  “hum,”  or  earth’s  background  free  oscillation  (Nishida  et  al.,  2009).  Furthermore,  we  extend 
the  useful  period  from  400  s  (Nishida  et  al.,  2009)  up  to  600  s!  Notice  the  remarkably  good  symmetry  of  very-long 
period  (>100  s)  EGFs  for  the  positive  and  negative  time  lags.  These  long-period  waves  help  extend  surface  wave 
constraints  to  all  depths  of  the  upper  mantle  and  uppermost  lower  mantle. 


Figure  5.  Example  of  EGFs  derived  ambient  noise  (Earth  “hum”)  with  various  wave  periods  (upper  left: 

50-600  s;  upper  right:  100-600  s;  lower  left:  200-600  s;  and  lower  right:  300-600  s).  The  lines  in  the 
accompanying  maps  show  the  paths  with  high  signal-to-noise  arrivals  from  the  virtual  source 
(station  HYB)  to  other  receivers  in  the  eastern  hemisphere.  The  red  dots  mark  a  wave  propagating 
at  a  velocity  of  3.9  km/s. 

We  measure  Rayleigh  wave  phase  delays  by  cross-correlating  observed  and  synthetic  Green’s  functions.  The 
waveforms  are  band-pass  filtered  in  six  period  ranges  (50-100,  75-150,  100-200,  150-300,  200-400,  and  300-600  s). 
A  total  of  ~1 1,000  phase  delay  measurements  meet  our  data  selection  criteria  (waveform  correlation  coefficient 
greater  than  0.85  and  signal-to-noise  ratio  greater  than  7  for  continental  station  pairs  and  5  for  those  with  at  least  one 
island  station).  Most  of  the  propagation  paths  with  useful  waves  are  in  the  northern  hemisphere  (Figure  6). 
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Figure  6.  Great  circle  paths  of  useful  Rayleigh  wave  delays  from  EGFs  at  200-400  (top),  100-200  (middle), 
and  50-100  s  (bottom)  periods.  Triangles  mark  the  locations  of  seismic  stations.  Eurasia,  northern 
Africa  and  Australia  are  well  sampled  at  all  periods  used.  The  Indian  and  Atlantic  oceans  are 
poorly  sampled,  especially  at  relatively  short  (50-100  s)  periods. 


Vp.  IwiiDolnl  vrt  nl  ZS3.59  ton  depth  Vi,  Irodrontnl  nirf  at  5U  ton  depth 


Figure  7.  Finite-frequency  sensitivity  of  long-period  EGFs  between  BJT  (Beijing)  and  LSA  (Lhasa)  to  Vp 

(left)  and  Vs  (right)  perturbations.  The  Rayleigh  wave  period  is  150-300  s.  The  horizontal  slice  is  at 
-250  km  depth.  Notice  that  the  sensitivity  to  Vp  is  mostly  in  the  upper  100  km.  The  unit  of  the 
sensitivity  kernel  is  s/m3. 

Finite-Frequency  Sensitivity  Kernels 

The  cross-correlation  of  vertical-vertical  component  yields  predominantly  Rayleigh  waves  in  the  time  window  and 
frequency  of  our  interest.  Although  commonly  neglected  in  surface  wave  tomography,  P-wave  speed  affects 
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Rayleigh  waves  (Figure  7).  The  sensitivity  to  P-wave  speed  is  near  the  surface  (the  depth  of  sensitivity  increases 
with  period)  and  at  a  level  comparable  to  the  Vs  kernel  near  the  surface. 


Inversion  and  the  Preliminary  Velocity  Model 

We  represent  the  Rayleigh  wave  delays  as  a  joint  Vp  and  Vs  inverse  problem  linearized  relative  to  an  iteratively 
updated  3D  model, 

<5t=J  [Ka{m0,x)Ama+K^(m0,x)Am/3]dV  (1) 


where  8t  is  the  travel  time  delay  relative  to  the  synthetics  for  the  3D  model  m0,  Ama  and  Amp  perturbations  to 
P-wave  and  S-wave  speeds,  respectively,  Ka  and  the  sensitivity  kernels  to  P-wave  and  S-wave  speeds, 
respectively.  The  integration  is  for  the  3D  volume  of  the  model. 
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Figure  8.  (Top)  A  horizontal  slice  of  the  4th-iteration  Vs  model  at  236  km  depth  from  full-wave  ambient  noise 
tomography.  The  color  scale  is  velocity  perturbation  (%)  relative  to  the  average  model.  The  plate 
boundaries  are  shown  as  dotted  lines.  The  2000-m  elevation  contour  (brown)  approximately  outlines 
the  Tibetan  plateau  and  other  mountain  ranges.  (Bottom)  The  Vs  velocity  structure  of  CUB  2.0 
(Ritzwoller  et  al.,  2002)  for  comparison. 
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The  above  equation  is  discretized  into  inversion  cells.  In  the  first  model  iterations  for  the  eastern  hemisphere,  we  use 
2°x2°  blocks  as  our  inversion  grids.  The  thickness  of  the  blocks  changes  from  -14  km  near  the  surface  to  -80  km  at 
-1600  km  depth.  The  inversion  block  size  may  be  reduced  in  late  iterations  to  capture  sharp  velocity  changes.  The 
inverse  problem  is  solved  with  LSQR  (Paige  and  Saunders,  1982)  with  damping  and  smoothness  constraints.  The 
damping  and  smoothness  constraints  are  gradually  reduced  until  the  data  misfit  is  consistent  with  the  estimated  data 
uncertainty  (normalized  x2~l),  as  in  Montelli  et  al.  (2004). 


Figure  8  is  a  comparison  of  the  4th- iteration  Vs  model  from  this  study  and  that  of  Ritzwoller  et  al.  (2002),  which 
was  derived  from  group  and  phase  velocities  of  earthquake  surface  wave  arrivals.  There  are  notable  similarities  and 
differences  between  the  two  models.  For  example,  the  African  cratons  and  the  west  Australian  craton  exhibit  high 
velocity  anomalies  in  both  models.  The  Afar  and  Eastern  Africa  Rift  are  prominent  low  velocity  anomalies  in  both 
models.  One  major  difference  between  the  two  models  is  the  structure  beneath  the  Tibetan  plateau,  an  area  well 
sampled  by  ambient  noise  (Figure  6).  CUB  2.0  has  a  strong  high-velocity  anomaly  beneath  the  center  of  the  plateau 
and  low-velocity  anomalies  in  the  surrounding  areas.  In  contrast,  full-wave  ambient  noise  tomography  reveals 
high-velocity  anomalies  associated  with  the  underthrusting  Indian  lithosphere  in  the  south,  the  Yangtzi  craton  in  the 
east  and  the  Sino-Korean  craton  in  the  north.  Vs  is  relatively  low  beneath  the  central  and  northern  plateau,  an 
observation  supported  by  regional  body  wave  tomography  (Li  et  al.,  2008). 
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Vp  (depth  =  78  km) 
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Figure  9.  The  4th-iteration  Vs  (top)  and  Vp  (bottom)  velocity  structures  at  78  km  depth.  The  color  scales  are 
different  for  Vs  and  Vp. 


One  of  the  reasons  for  the  differences  between  the  two  models  in  Figure  8  lies  in  the  fact  that  phase  delays  are 
attributed  to  both  Vs  and  Vp  perturbations  in  full-wave  tomography  (Figure  9).  The  similarity  in  major  features 
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between  Vs  and  Vp  in  Figure  9  suggest  that  Rayleigh  waves  yield  useful  resolution  of  the  Vp  structure  in  the  crust 
and  shallow  mantle. 


The  very  long  period  EGFs  allow  us  to  extend  resolution  through  the  entire  upper  mantle  (Figure  10).  The  mantle 
transition  zone  beneath  the  Afar  and  East  African  Rift  has  a  strong  low  velocity,  suggesting  that  that  rifting  and 
hotspot  volcanism  have  a  source  rooted  at  least  in  the  mantle  transition  zone.  There  is  a  strong  high-velocity 
anomaly  immediately  south  of  Tibet.  This  is  interpreted  as  evidence  for  the  Tethyan  subducted  slabs  under  India 
(Van  der  Voo  et  al.,  1999).  Perhaps  the  most  intriguing  result  in  the  transition  zone  is  the  strong  low  velocity 
anomalies  north  of  the  subducted  Tethyan  slabs  (Hafkenscheid  et  al.,  2006).  One  possible  explanation  is  that  these 
anomalies  represent  elevated  water  content  released  from  the  subducted  Tethyan  slabs  (Nolet  and  Zielhuis,  1994). 
van  der  Lee  et  al.  (2008)  observed  a  similar  low-velocity  anomaly  in  the  upper  mantle  transition  zone  beneath  the 
U.S.  East  Coast  and  attributed  it  to  water  released  from  the  subducted  Farallon  lithosphere. 
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Figure  10.  Vs  structure  at  521  km  depth  from  full-wave  ambient  noise  tomography.  The  color  scale  is  velocity 
perturbation  (%)  relative  to  an  average  model. 


CONCLUSIONS 


We  collected  and  processed  up  to  21  years  continuous  seismic  data  from  184  broadband  seismic  stations  in  the 
eastern  hemisphere  (latitude  55°S-55°N;  longitude  30°W-157°E).  EGFs  derived  from  ambient  seismic  noise  show 
clear  Rayleigh  waves  over  a  wide  frequency  range.  We  measured  ~1 1,000  phase  delays  in  6  period  bands  ranging 
from  50  s  to  600  s.  The  very  long  period  Rayleigh  waves  extend  Vs  resolution  through  the  entire  upper  mantle.  The 
new  model  provides  constraints  on  several  geological  processes  in  the  region,  such  as  the  origin  of  the  East  African 
Rift,  the  extent  of  underthrusting  of  the  Indian  lithosphere  beneath  Tibet,  and  recycling  of  water  by  subducted  slabs. 
Each  of  these  geological  processes  warrants  detailed  studies  that  will  be  accomplished  by  the  integration  of  ambient 
noise  and  earthquake  data  and  the  development  of  a  hierarchical,  multi-grid  (resolution)  model. 
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